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Abstraction  in  Numerical  Methods 


Matthew  Halfant  and  Gerald  Jay  Sussman 


Abstract 


We  illustrate  how  the  liberal  use  of  high-order  procedural  abstrac¬ 
tions  and  infinite  streams  helps  us  to  express  some  of  the  vocabulary 
and  methods  of  numerical  analysis.  We  develop  a  software  toolbox  en¬ 
capsulating  the  technique  of  Richardson  extrapolation,  and  we  apply 
these  tools  to  the  problems  of  numerical  integration  and  differentia¬ 
tion.  By  separating  the  idea  of  Richardson  extrapolation  from  its  use 
in  particular  circumstances  we  indicate  how  numerical  programs  can 
be  written  that  exhibit  the  structure  of  the  ideas  from  which  they  are 
formed . 


Keywords:  Scheme,  Abstraction,  Programming  Methodology,  Richard¬ 
son  Extrapolation,  LISP. 
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Abstract 

We  illustrate  how  the  liberal  use  of  high-order  procedural  abstrac¬ 
tions  and  infinite  streams  helps  us  to  express  some  of  the  vocabulary 
and  methods  of  numerical  analysis.  We  develop  a  software  toolbox  en¬ 
capsulating  the  technique  of  Richardson  extrapolation,  and  we  apply 
these  tools  to  the  problems  of  numerical  integration  and  differentia¬ 
tion.  By  separating  the  idea  of  Richardson  extrapolation  from  its  use 
in  particular  circumstances  we  indicate  how  numerical  programs  can 
be  written  that  exhibit  the  structure  of  the  ideas  from  which  they  are 
formed. 

A  numerical  analyst  uses  powerful  ideas  such  as  Richardson  extrapolation 
for  organizing  programs,  but  numerical  programs  rarely  exhibit  the  structure 
implied  by  the  abstractions  used  in  their  design.  It  is  traditional  practice  in 
the  domain  of  numerical  methods  for  each  program  to  be  hand  crafted,  in 
detail,  for  the  particular  application,  rather  than  to  be  constructed  by  mixing 
and  matching  from  a  set  of  interchangable  parts.  Such  numerical  programs 
are  often  difficult  to  write  and  even  more  difficult  to  read. 

In  this  paper  we  illustrate  how  the  liberal  use  of  high-order  procedural 
abstractions  and  infinite  streams  helps  us  to  express  some  of  the  vocabulary 
and  methods  of  numerical  analysis.  We  develop  a  software  toolbox  encapsu¬ 
lating  the  technique  of  Richardson  extrapolation,  and  we  apply  these  tools  to 
the  problems  of  numerical  integration  and  differentiation.  By  separating  the 
idea  of  Richardson  extrapolation  from  its  use  in  particular  circumstances  we 
indicate  how  numerical  programs  can  be  written  that  exhibit  the  structure 
of  the  ideas  from  which  they  are  formed. 


Figure  1:  |ABj  =  5n,  |AC|  =  |CB|  =  52n 


A  first  example:  Archimedean  computation  of  x 

We  begin  with  a  playful  example:  approximating  the  value  of  x  by  the 
method  of  Archimedes.  Let  5„  be  the  length  of  one  side  of  a  regular  n- 
sided  polygon  inscribed  in  a  unit  circle.  As  n  approaches  infinity,  the  semi¬ 
perimeter  Pn  =  nSn/ 2  approaches  x.  Applying  the  Pythagorean  Theorem  to 
right  triangles  ACD  and  ADO  from  figure  1  we  derive  the  relation: 

SJn  =  J2-J4-  S’. 

Equivalently, 

o  ^ n 

*-,2rj  /  =  ? 

the  latter  form  being  preferred  because  it  avoids  the  subtraction  of  nearly 
equal  quantities  as  Sn  — >  0.  In  Scheme  [4]  (the  dialect  of  Lisp  we  use)  we 
can  write  the  transformation  from  Sn  to  S^n  as: 

(define  (ref ine-by-doubling  s)  ;  s  is  a  side 

(/  s  (sqrt  (+  2  (sqrt  (-  4  (*  s  s))))))) 


Starting  with  a  square  (54  is  simply  we  want  to  form  the  sequence 
of  side  lengths  54, 5s,  5 16,  ■  •  ••  Such  sequences  are  naturally  represented  with 
streams — effectively  infinite  lists  whose  terms  are  evaluated  only  on  demand 


f 

r 


(for  a  discussion  of  streams  see  [1]).  We  use  a  stream  generator  to  produce 
the  orbit  of  a  starting  value  under  repeated  application  of  a  transformation 
next: 

(define  (stream-of-iterates  next  value) 

(cons-stream  value 

(stream-of -iterates  next  (next  value)))) 

Now  we  can  define  the  stream  of  side  lengths: 

(define  side-lengths 

(stream-of-iterates  refine-by-doubling  (sqrt  2))) 
and  the  corresponding  stream  of  numbers  of  sides: 

(define  side-numbers 

(stream-of-iterates  (lambda  (n)  (*  2  n))  4)) 

Combining  these  termwise  using  map-streams  lets  us  form  the  sequence  of 
semi-perimeters  P4,P8,.Pi6,  •  •  *,  whose  limit  is  it: 

(define  (semi-perimeter  length-of-side  number-of-sides) 

(*  (/  number-of-sides  2)  length-of-side)) 

(define  archimedean-pi-sequence 

(map-streams  semi-perimeter  side-lengths  side-numbers)) 

We  can  look  at  the  results: 

(print-stream  archimedean-pi-sequence) 

==>  2.82842712474619 

3.06146745892072 
3.12144515225805 
3.13654849054594 
3.14033115695475 
3.14127725093277  ;  term  6 

3.14159265358862 
3.14159265358950 
3.14159265358972 
3.14159265358977 


;  term  20 


3.14159265358979 

3.14159265358979 

3.14159265358979 

As  expected,  the  sequence  converges  to  7r,  but  it  takes  24  terms  to  reach 
the  full  machine  precision.  Imagine  poor  Archimedes  doing  the  arithmetic 
by  hand:  square  roots  without  even  the  benefit  of  our  place  value  system! 
He  would  be  interested  in  knowing  that  full  precision  can  be  reached  on  the 
fifth  term,  by  forming  linear  combinations  of  the  early  terms  that  allow  the 
limit  to  be  seized  by  extrapolation. 

To  understand  how  this  is  done,  rewrite  the  expression  for  the  semi¬ 
perimeter  Pn  by  using  the  Taylor  series  for  Sn  =  2  sin  (x /n)  (see  figure  1): 

Pn  =  (n/2  )Sn 

=  (n/2)(2sin  (7r/n)) 

=  7 r  +  A/n 2  +  B/ii'  +  •  •  • 

where  A  and  B  are  constants  whose  values  aren’t  important  here.  As  n  gets 
larger,  it  is  the  A/n 2  term  that  dominates  the  truncation  error — the  differ¬ 
ence  between  7r  and  Pn  for  finite  n.  Whenever  we  double  n,  this  principal 
component  of  the  truncation  error  is  reduced  by  a  factor  of  4;  that  knowl¬ 
edge  allows  us  to  combine  successive  terms  of  the  sequence  •  •  •  Pn,P2n,  to 
eliminate  entirely  the  effect  of  this  quadratic  term. 

Specifically:  multiply  P2n  by  4,  so  that  its  l/n2  term  matches  that  of  Pn ; 
then  subtract  Pn  to  eliminate  this  term;  after  some  slight  rearrangement  we 
get 

4  P2n  ~Pu  B 

-  —  TT - -  -f  •  •  • 

3  4n4 

The  accelerated  sequence  P un  is  defined  by  the  expression  on  the  left;  it 
approximates  x  with  a  truncation  error  that  goes  to  0  like  l/n4,  which  is 
much  faster  than  before. 

We  could  revise  our  program  to  compute  the  sequence  PA,  P{6,  •  ■  •,  and 
demonstrate  that  its  convergence  is  more  rapid  than  that  of  P4,  P8,  Pie,  •  • 
However,  this  acceleration  scheme  is  quite  general,  and  we  prefer  to  develop 
it  away  from  the  present  context.  Afterwards,  we  can  apply  it  to  this  and 
other  pursuits. 


The  method  of  Richardson  extrapolation 

Instead  of  phrasing  our  argument  in  terms  of  a  parameter  n  that  gets 
larger  through  successive  doubling,  we’ll  follow  the  convention  of  using  a 
positive  quantity  h  that  approaches  0  through  successive  divisions  by  2.  In 
the  example  above,  we  need  only  make  the  identification  h  =  I/n.  Sometimes 
n,  sometimes  h,  will  be  the  more  natural,  and  we  can  formulate  the  problem 
however  we  prefer.  Typically  h  will  arise  as  a  step  size. 

Now  imagine  that  we  seek  the  limit,  as  h  — *  0,  of  some  function  R(h), 
and  that  we  pursue  it  by  constructing  the  sequence  R(h),  R(h/2),  R(h/ 4),  •  •  •. 
We  suppose  that  each  expression  R{h)  represents  the  limiting  value  A  with 
a  truncation  error  that  is  analytic  at  h  =  0: 

R(h)  =  /l  +  Bhp'  +  Chp 2  +  Dhn  +  •  •  • 

The  exponents  Pi,P2,P3,- '  •  may  well  be  just  the  natural  numbers  1,2,3,-  •  •; 
but  they  might  represent  a  subsequence — for  example,  the  even  numbers, 
as  they  did  in  the  Archimedean  example  above.  We  assume  this  sequence 
is  known.  On  the  other  hand,  we  need  not  know  in  advance  the  values  of 
A,  B,  C,  •  •  •;  indeed,  our  whole  purpose  is  to  determine  A. 

Knowing  the  truncation  error  as  a  power  series  in  h  allows  us  to  eliminate 
the  effect  of  the  dominant  term:  we  do  this  by  subtracting  the  appropriate 
multiple  of  R(h/ 2)  from  R(h): 


or 


R(h)  -2PlR(h/2)  =  (1  -  2 p')A  +  ClhP2  +  Dxhp>  + 
2p'R{h/2)  -  R(h) 


R\h)  = 


2  pi  -  1 
=  A  +  C2hP2  +  D-jhP3  + 


We  view  this  as  a  transformation  applied  to  (adjacent  terms  of)  the  se¬ 
quence 

R(h),  R(h/2),  R(h/4),  ■  ■  ■ 
to  produce  a  new  sequence 


V 


Figure  2:  The  Richardson  Tableau 


whose  truncation  error  is  now  dominated  by  the  term  containing  h p2.  That 
is,  the  R'  sequence  converges  faster  than  the  original. 

Now  we  come  to  Richardson’s  great  idea:  since  R'  has  a  truncation  error 
dominated  by  C2/11'2 ,  we  can  apply  the  same  idea  again: 


R"(h)  = 


2P2 R'(n/2)  -  R’(h) 
-  1 


A  +  D^hV3  +  •  ■  • 


yielding  a  sequence  R"(h ),  R"(h/ 2),  R"(h/4),  ■  ■  •  whose  truncation  error  now 
has  the  dominant  term  £>3/iP3.  And  so  on.  Given  the  sequence  pi,P2,P3,  •  •  •, 
one  can  form  a  tableau  (see  figure  2)  in  which  the  original  sequence  appears 
as  the  vertical  s  column  at  the  left;  to  the  right  is  the  derived  t  column;  the 
u  column  is  derived  from  t  as  t  is  from  s,  and  so  on.  All  columns  converge  to 
the  same  limit,  /f(0+),  but  each  converges  faster  than  its  predecessor.  Thus, 
s  converges  with  an  error  term  0(hPl) — that  is,  of  order  hPi ;  t  converges 
with  error  term  0{hP2)\  u  has  error  term  0(hP3),  and  so  on. 

By  the  way,  if  the  sequence  Pi,P2,P3, •••  is  not  known  in  advance,  one 
can  take  the  conservative  approach  of  assuming  it  to  be  the  natural  number 
sequence  1,  2, 3,  •  •  •;  this  leads  at  worst  to  the  inefficiency  of  creating  adjacent 
columns  with  the  same  error  order.  Alternatively,  the  appropriate  exponent 
p  for  a  given  column  can  be  inferred  numerically  from  the  early  terms  of  that 
column;  we  have  done  that  in  our  library,  but  do  not  pursue  it  here. 


The  Richardson  Toolbox 

Having  sketched  the  basic  ideas,  we  develop  Richardson  extrapolation  as 
a  set  of  tools  that  can  be  applied  in  diverse  contexts.  First,  if  we  don’t  already 
have  the  sequence  R(h),  R(h/2),  R(h/4),  ■  ■  ■  we  need  to  be  able  to  make  it, 


given  R  and  h.  The  name  make-zeno-sequence  derives  from  the  suggestive 
connection  to  successive  halving  in  the  statement  of  Zeno’s  paradox. 


(define  (make-zeno-sequence  R  h) 

(cons-stream  (R  h)  (make-zeno-sequence  R  (/  h  2)))) 

The  basic  operation  of  accelerating  the  sequence  requires  that  we  know 
the  order  of  the  dominant  error  term.  Accelerate-zeno-sequence  takes 
this  order  as  an  explicit  argument. 

(define  (accelerate-zeno-sequence  seq  p) 

(let*  ( (2~p  (expt  2  p))  (2"p-l  (-  2~p  1))) 

(map-streams  (lambda  (Rh  Rh/2) 

(/  (-  (*  2“p  Rh/2)  Rh) 

2-p-l)) 

seq 

(tail  seq)))) 

We  can  get  our  hands  on  the  full  tableau.  In  this  case  we  iterate  the 
application  of  accelerate-zeno-sequence  to  make  an  infinite  sequence  of 
accelerated  sequences.  In  forming  the  full  tableau,  we  make  the  simplifying 
assumption  that  the  exponent  sequence  is  arithmetic: 


{Pu  Pi,  7>3,  •••}  =  il>*  P  +  <Z,  P  i-  2<z,  p  +  3q ,•••} 

Hence  it  can  be  specified  by  the  initial  order  p  and  an  increment  q.  In  typical 
cases,  q  is  1  or  2. 

The  procedure  make-zeno-tableau  accepts  as  arguments  the  original 
Zeno  sequence  along  with  the  characterizing  order  p  and  increment  q ;  it 
returns  the  sequence  of  accelerated  sequences. 

(define  (make-zeno-tableau  seq  p  q) 

(define  (sequences  seq  order) 

(cons-stream  seq 

(sequences  (accelerate-zeno-sequence  seq  order) 

(+  order  q)))) 

(sequences  seq  p)) 
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Finally,  the  procedure  f irst-terms-of-zeno-tableau  produces  a  se¬ 
quence  si,  tl,  ul,  vl,  wl,  ...  of  t  he  first  terms  taken  from  each  of  t  he 
accelerated  sequences.  T  he  sequence  ol  first  terms  converges  at  a  remarkable 
rate:  not  only  are  the  first  n  terms  of  the  original  truncation-ei  ror  series 
removed,  but  the  remainder  is  effectively  divided  by  2”  '^2  (see  appendix  for 
details). 

(define  (f irst-terms-of-zeno-tableau  tableau) 

(map-streams  head  tableau)) 

(define  (richardson-sequence  seq  p  q) 

(f irst-terms-of-zeno-tableau  (make-zeno-tableau  seq  p  q))) 


Archimedes  revisited 

Before  proceding  any  further,  let's  see  how  well  this  does  on  our  original 
example.  We  apply  richardson-sequence  to  the  archimedean-pi-sequence 
previously  computed  and  examine  the  result: 

(print-stream  (richardson-sequence  archimedean-pi-sequence  2  2)) 

==>  2.82842712474619 

3.13914757031223 
3.14159039312994 
3.14159265328605 

3.14159265358979 

3.14159265358979 

3.14159265358979 


As  indicated  earlier,  full  precision  is  reached  on  the  5th  term  (although 
we  need  to  compute  the  6th  to  know  that  we’ve  reached  it,  and  may  want 
to  compute  the  7th  just  to  be  sure).  What  remains  now  is  to  establish  some 
algorithmic  means  for  assert.aining  when  a  limit  has  been  reached. 

Completing  and  applying  the  Richardson  toolbox 

Now  that  we  have  the  idea  under  control,  we  must  fdl  in  our  Richardson 
toolbox  to  allow  its  application  in  a  variety  of  situations.  We  need  ways  of 
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extracting  our  best  estimate  of  the  limit  from  a  sequence.  One  simple  crite¬ 
rion  that  may  be  used  is  this:  We  declare  convergence  when  two  consecutive 
terms  are  sufficiently  close. 

Alas,  the  notion  of  sufficient  closeness  is  slightly  sticky:  relative  accuracy 
is  what’s  generally  wanted,  but  that  fails  in  the  case  of  a  sequence  with  limit 
0.  One  way  around  the  difficulty  is  to  use  the  metric 

\hi-h2\ 

( |At !  +  |  /t  2 1 )  /  2  +  1 

to  measure  the  distance  between  two  numbers  h\  and  /12: 

(define  (close-enuf?  hi  h2  tolerance) 

(<=  (abs  (-  hi  h2) ) 

(*  .5 

tolerance 

(+  (abs  hi)  (abs  h2)  2)))) 

This  criterion  amounts  to  relative  closeness  when  the  numbers  to  be  com¬ 
pared  are  large,  but  makes  a  graceful  transition  to  absolute  closeness  when 
the  numbers  are  much  smaller  (in  magnitude)  than  1. 

Using  this  or  any  similar  predicate,  we  construct  our  limit  detector: 

(define  (stream-limit  s  tolerance) 

(let  loop  ((s  s)) 

(let*  ((hi  (head  s))  (t  (tail  s))  (h2  (head  t))) 

(if  (close-enuf?  hi  h2  tolerance) 
h2 

(loop  t))))) 

A  more  cautious  version  of  the  limit  detector  would  require  close  agree¬ 
ment  for  three  or  more  successive  terms  (we’ve  been  bitten  ourselves  by 
accidental  equality  of  the  first  two  terms  of  a  sequence).  Actually,  there’s 
another  modification  we’ll  be  forced  to  make  very  shortly:  we’ll  need  an 
optional  final  argument  m  that  forbids  stream-limit  from  examining  more 
than  the  first  m  terms  of  the  sequence  before  returning  an  answer. 

Given  stream-limit,  the  following  combination  proves  useful  for  finding 
the  Richardson  limit  of  a  function.  The  arguments  ord  and  inc  are  our 
previous  p  and  </. 


(define  (richardson-limit  f  start-h  ord  inc  tolerance) 

(stream-limit 

(richardson-sequence  (make-zeno-sequence  f  start-h) 

ord 

inc) 

tolerance) ) 

We  are  ready  now  to  apply  our  tools  to  a  significant  example. 

Numerical  computation  of  derivatives 

The  following  higher-order  procedure  takes  a  procedure  that  computes  a 
numerical  function  and  returns  a  procedure  that  calculates  an  approximation 
to  the  derivative  of  that  function: 

(define  (make-derivative-function  f) 

(lambda  (x) 

(let  ( (h  .00001)) 

(/  (-  (f  (+  x  h))  (f  (-  x  h))) 

2  h)))) 

Notice  the  ad  hoc  definition  of  h.  We  are  walking  the  line  between  trunca¬ 
tion  error  (not  having  h  small  enough  for  the  difference  quotient  to  adequately 
approximate  the  derivative)  and  roundoff  error  (having  h  so  small  that  the 
subtraction  of  nearly  equal  quantities  loses  all  accuracy  in  the  answer).  The 
optimal  h  depends  both  on  x  and  on  the  number  of  digits  carried  by  the  ma¬ 
chine,  but  even  with  this  h  we’ll  generally  lose  about  a  third  of  our  significant 
digits  (we’d  lose  half  of  our  digits  had  wc  used  the  forward,  rather  than  the 
centered  difference  quotient).  Of  course  we’re  hoping  that  Richardson  will 
allow  us  to  do  better. 

It’s  instructive  to  experiment  with  letting  h  go  to  0.  Given  f ,  x,  and  h,  we 
produce  a  stream  of  difference  quotients  in  which  h  is  successively  reduced 
by  a  factor  of  2. 

(define  (dif f-quot-stream  f  x  h) 

(cons-stream  (/  (-  (f  (+  x  h))  (f  (-  x  h)))  2  h) 

(diff-quot-stream  f  x  (/  h  2)))) 

We  apply  this  to  the  estimation  of  the  derivative  of  the  square  root  at  1 
(exact  answer  is  0.5). 


(print-stream  (diff-quot-stream  sqrt  1  .1)) 


==>  0.500627750598189 

0.500156421150633 
0.500039073185090 
0.500009766292631 
0.500002441447984 
0.500000610354192 
0.500000152587994 
0.500000038146951 
0.500000009536734 
0.500000002384411 
0.500000000595833 
0.500000000148475 
0.500000000038199 

0.500000000010914  ;  14th  term 

0.500000000010914 

0.499999999974534 

0.499999999992724 

0.500000000029104 

0.499999999883585 

0.500183105468750  ;  40th  term 

0.500488281250000 

0.499267578125000 

0.498046875000000 

0.502929687500000 

0.507812500000000 

0.488281250000000 

0.468750000000000 

0.546875000000000 

0.625000000000000 

0.312500000000000  ;  50th  term 

0. 

0. 

0. 


We  observe  that  the  error  diminishes  steadily  until  the  14th  term  is 
reached;  after  this,  the  error  builds  back  up  in  a  somewhat  erratic  man- 
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ner  until,  after  the  50th  term,  we  are  left  with  a  steady  parade  of  zeros. 
This  problem  results  from  the  subtraction  of  nearly  equal  quantities  in  the 
numerator  of  the  difference  quotient:  we  lose  more  and  more  significant  fig¬ 
ures  until  h  becomes  so  small  that  x  +  h  and  x  —  h  are  equal  to  full  working 
precision,  after  which  only  0  quotients  can  be  returned. 

Hence  we  are  in  a  race  between  truncation  error,  which  starts  out  large 
and  gets  smaller,  and  roundoff  error,  which  starts  small  and  gets  larger. 
Richardson  helps  the  situation  by  creating  new  sequences  in  which  the  trun¬ 
cation  error  dimishes  more  rapidly,  which  is  just  what  we  need.  To  be  more 
precise,  we  need  to  look  at  how  the  roundoff  error  works  in  this  example. 

Any  real  number  x,  represented  in  the  machine,  is  rounded  to  a  value 
2(1  +  e),  where  e  is  effectively  a  random  variable  whose  absolute  value  is  on 
the  order  of  the  machine  epsilon ,  e:  that  smallest  positive  number  for  which 
1.0  and  1.0  +  e  can  be  distinguished.  For  IEEE  double  precision  (as  used,  for 
example,  by  the  8087  numeric  coprocessor),  e  =  2-53  =  1.11  x  10~16.  Now 
if  h  is  small,  both  f(x  +  h)  and  f(x  —  h)  have  machine  representations  in 
error  by  around  f(x)r,  their  difference  suffers  an  absolute  error  of  this  same 
order.  Since  the  difference  f(x  +  h)  —  f(x  —  h)  should  equal  around  f'(x)2h, 
the  relative  error  is  of  the  order 

f  /(*)  I 

The  relative  error  of  the  difference  quotient  is  essentially  the  same  as  that  of 
its  numerator,  the  denominator  being  just  2h  which  is  known  to  full  precision. 

From  the  above  expression,  we  see  that  the  relative  error  due  to  roundoff 
basically  doubles  each  time  h  is  halved — a  result  that  is  easy  to  see  directly 
in  terms  of  the  binary  representation  of  x  +  h:  dividing  h  by  2  shifts  the 
binary  representation  of  h  one  position  to  the  right;  but  the  presence  of  x 
nails  down  the  high  order  bits  of  x  +  /;,  whence  the  low  order  bits  of  h  fall 
off  the  end,  one  per  iteration. 

Suppose  we  want  to  compute  the  derivative  of  the  square  root  at  1  with 
a  relative  error  of  at  most  10-13,  and  starting  with  h  =  0.1.  We  need  to 
estimate  the  initial  relative  roundoff  error;  the  preceding  formula  must  be 
modified  slightly  for  this  purpose.  First,  the  denominator  is  actually  f(x  + 
h)  —  f(x  —  h),  which  is  what  we  must  use  (it  was  written  above  as  2 hf(x) 
only  to  show  the  trend  as  h  gets  small).  Second,  we  want  to  ensure  that 


the  predicted  relative  roundoff  error  is  at  least  a  positive  multiple  of  the 
machine  epsilon;  hence  we  take  the  next-highest-integer  of  the  absolute- value 
subexpression: 

vT 

vn~  n/09 

Thus,  the  initial  relative  roundoff  error  is  10  roundoff  units ,  or  10c  =  1.1  x 
10~15.  Since  the  roundoff  error  roughly  doubles  at  each  iteration,  we  ask: 
How  many  times  can  1.1  x  10-15  be  doubled  before  reaching  10~13? 


)  =  10 


1  +  floor( 


1.1  x  10"152n  <  10~13 


so 


log  (10_13/1.1  x  10-15) 
log  2 


=  6.5 


Hence  if  we  restrict  ourselves  to  at  most  6  terms  past  the  first  (for  a  total  of 
7  terms),  we  can  be  reasonably  sure  our  data  is  uncontaminated  by  noise  at 
the  level  of  interest.  This  makes  accelerated  convergence  really  crucial:  we 
have  to  reach  our  limit  quickly  or  not  at  all. 

Here’s  a  modified  version  of  stream-limit  that  accepts  an  optional  nnal 
parameter  m,  designating  the  maximum  number  of  stream  terms  to  examine 
in  the  search  for  the  limit.  If  m  is  reached  without  convergence,  we  just 
return  the  final  term  as  a  best  guess;  a  more  professional  approach  would 
be  to  return  some  kind  of  an  error  code,  along  with  the  best  guess  and  an 
estimate  of  its  truncation  error. 


(define  (stream-limit  s  tolerance  .  opts)  ;opts  =  optional  args 
(let  ((M  (if  (null?  opts)  ’nomax  (car  opts)))) 

(let  loop  ((s  s)  (count  2)) 

(let*  ((hi  (head  s))  (t  (tail  s))  (h2  (head  t))) 

(if  (close-enuf?  hi  h2  tolerance) 
h2 

(if  (and  (number?  M)  (>=  count  M)) 
h2 

(loop  t  (+  count  1)))))))) 


The  revised  version  of  richardson-limit  is  simply: 


(define  (richardson-limit  f  start-h  ord  inc  tolerance  .  opts) 
(stream-limit 

(richardson-sequence  (make-zeno-sequence  f  start-h) 

ord 

inc) 

tolerance 

(if  (null?  opts)  ’nomax  (car  opts)))) 

We  can  now  define  our  derivative  estimator  in  the  following  natural  way.  (In 
practice,  the  routine  as  shown  admits  of  several  pitfalls:  h  becomes  0  if  x  is; 
delta  might  end  up  as  0  by  chance;  and  possibly  some  other  bad  things  we 
haven’t  thought  of.  It  will  serve  for  purposes  of  illustration  here.) 

(define  rderiv 

(lambda  (f  tolerance) 

(lambda  (x) 

(let*  ((h  (*  0.1  (abs  x))) 

(delta  (-  (f  (+  x  h))  (f  (-  x. h)))) 

(roundoff  (*  *machine-epsilon* 

(+  1  (floor  (abs  (/  (f  x)  delta)))))) 

(n  (floor  (/  (log  (/  tolerance  roundoff)) 

(log  2))))) 

(richardson-limit  (lambda  (dx) 

(/  (-  (f  (+  x  dx)) 

(f  (-  x  dx))) 

(*  2  dx))) 
h 
2 
2 

tolerance 
(+  n  1)))))) 

Notice  that  the  ord  and  inc  arguments  are  both  2:  the  truncation  error 
involves  only  even  powers  of  h.  Had  we  used  the  forward  difference  quotient, 
(/(x  +  h)  —  f(x))/h,  then  all  powers  of  h  would  arise,  and  ord  and  inc  would 
both  be  1.  These  results  follow  easily  from  the  Taylor  expansion  of  /. 

Applied  to  the  square  root  example,  we  find: 

((rderiv  sqrt  le-13)  1)  ==>  0.500000000000016 


which  shows  a  relative  error  of  0.32  x  10-13.  Further  testing  shows  that  the 
relative  error  of  10-13  is  generally  met. 

We  pass  now  to  another  significant  application,  in  which  roundoff  error 
is  happily  not  an  issue. 

Numerical  integration  by  Romberg’s  method 

Given  a  function  /  that  behaves  nicely  (i.e.,  has  two  continuous  deriva¬ 
tives)  on  a  finite  interval  [a,  6],  we  seek  numerical  approximations  to  the 
definite  integral  of  /  from  a  to  b.  The  plan  of  attack  is  to  divide  [a,  6]  into 
some  number,  n,  of  subintervals  each  of  length  h  =  (6  —  a)/n.  We  apply 
the  trapezoidal  rule  to  compute  an  approximating  sum  Sn;  we  then  form  a 
sequence  of  approximations  by  repeatedly  doubling  n  (equivalently,  halving 
h)  and  use  Richardson  extrapolation  on  the  result. 

We  get  the  ball  rolling  with  the  following  procedure.  It  takes  /,  a ,  and  6, 
and  returns  a  procedure  that,  given  n,  computes  S„ : 

(define  (trapezoid  fab) 

... .  (lambda  (n) 

(o"'  (let  ((h  (/  (-  b  a)  n))) 

(let  loop  ((i  1)  (sum  (/  (+  (f  a)  (f  b))  2))) 

(let  ((x  (+  a  (*  i  h)))) 

(if  «  i  n) 

(loop  (+  i  1)  (+  sum  (f  x))) 

(*  sura  h))))))) 

We  use  this  to  estimate  tt : 

(define  (f  x)  (/  4  (+  1  (*  x  x)))) 

(define  pi-estimator  (trapezoid  f  0  1)) 

(pi-estimator  10)  ==>  3.13992598890716 

(pi-estimator  10000)  ==>  3.14159265192314 

It  is  shown  in  standard  texts  (for  example  [2]  or  [3])  that,  for  /  in  C2[a,6] 
as  we’ve  assumed,  the  truncation  error  involves  only  even  powers  of  h.  Hence 
we  proceed: 

(define  (pi-estimator-sequence  n) 

(cons-stream  (pi-estimator  n) 

(pi-estimator-sequence  (*  2  n)))) 


(print-stream 

(richardson-sequence 

(pi-estimator-sequence  10)  22)) 

==>  3.13992598890716 

3.14159265296979 
3.14159265362079 
3.14159265358979 

The  convergence  rate  is  very  encouraging — we  get  full  machine  accuracy 
in  only  80  evaluations  of  the  original  function  /—  but  there  is  considerable 
redundant  computation  here.  Every  time  vve  double  the  number  of  points 
we  reevaluate  the  integrand  at  the  old  grid  points;  this  is  a  lamentable  in¬ 
efficiency  in  cases  where  /  is  expensive  to  compute.  Romberg’s  method  of 
quadrature,  to  which  we  now  proceed,  is  a  variation  of  the  above  that  avoids 
unnecessary  recomputation  of  /. 

We  begin  with  a  utility  procedure  that  computes  a  sum  of  terms  f(i), 
where  i  takes  unit  steps  from  a  to  (not-greater-than)  b : 

(define  (sigma  fab) 

(let  loop  ((sum  0)  (x  a)) 

(if  (>  x  b) 
sum 

(loop  (+  sum  (f  x))  (+  x  1))))) 

We  examine  how  the  sum  S2n  is  related  to  Sn.  In  the  former,  we  employ 
a  partition  of  [a,b]  into  2 n  equal  parts,  having  grid  points  at  Xi  =  a  +  ih, 
with  h  =  (b  —  a)/2n,  and  i  going  from  0  to  2 n.  The  even-index  terms  of 
this  point  sequence  comprise  the  entire  point  sequence  for  Sn;  in  computing 
S2n:  we  need  only  evaluate  /  at  the  odd  indices,  the  others  being  already 
incorporated  into  Sn: 

$2n  ~  ~Sn  +  h  f(x2i-\ ) 

1  «=1 

This  recursion  is  the  basis  of  the  following  procedure,  which  generates  the 
sequence  S\ ,  52 ,  S4 ,  •  ■  • : 


(define  (trapezoid-sums  fab) 

(define  (next-S  S  n) 

(let*  ((h  (/  (-  b  a)  2  n)) 

(fx  (lambda(i)  (f  (+  a  (*  (+  i  i  -1)  h) ) ) ) ) ) 

(+  (/  S  2)  (*  h  (sigma  fx  1  n))))) 

(define  (S-and-n-stream  S  n) 

(cons-stream  (list  S  n) 

(S-and-n-stream  (next-S  S  n)  (*  n  2)))) 

(let*  ( (h  (-  b  a)) 

(S  (*  (/  h  2)  (+  (f  a)  (f  b))))) 

(map-stream  car  (S-and-n-stream  S  1)))) 

We  arrive  at  the  more  economical  version  of  our  previous  method: 

(define  (romberg  fab  tolerance) 

(stream-limit 

(richardson-sequence  (trapezoid-sums  fab) 

2 

2) 

(0  tolerance)) 


Conclusion 

We  have  shown  how  a  classical  numerical  analysis  method,  Richardson 
extrapolation,  can  be  formulated  as  a  package  of  procedures  that  can  be  used 
as  interchangable  components  in  the  construction  of  traditional  applications 
such  as  the  estimation  of  a  derivative  and  Romberg  quadrature.  Such  a 
formulation  is  valuable  in  that  it  separates  out  the  ideas  into  several  inde¬ 
pendent  pieces,  allowing  one  to  mix  and  match  combinations  of  components 
in  a  flexible  way  to  facilitate  attacking  new  problems.  Clever  ideas,  such  as 
Richardson  extrapolation,  need  be  coded  and  debugged  only  once,  in  a  con¬ 
text  independent  of  the  particular  application,  thus  enhancing  the  reliability 
of  software  built  in  this  way.  Roylance  [5]  has  similar  goals.  He  constructs 
high-performance  implementations  of  special  functions,  abstracting  out  re¬ 
current  themes  such  as  Chebyshev  economization. 

The  decompositions  we  displayed  have  used  powerful  abstraction  mech¬ 
anisms  built  on  high-order  procedures,  interconnected  with  interfaces  orga¬ 
nized  around  streams.  The  programs  were  implemented  functionally — there 
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were  no  assignments  or  other  side-effects  in  any  of  our  example  programs. 
Because  functional  programs  have  no  side  effects  they  have  no  required  order 
of  execution.  This  makes  it  exceptionally  easy  to  execute  them  in  parallel. 

A  program  is  a  communication,  not  just  between  programmers  and  com¬ 
puters,  but  also  between  programmers  and  human  readers  of  the  program; 
quite  often,  between  the  programmer  and  him/herself.  A  program  describes, 
more  or  less  clearly,  an  idea  for  how  to  obtain  some  desired  results.  One 
power  of  programs  is  that  they  allow  one  to  make  the  knowledge  of  methods 
explicit,  so  that  methods  can  be  studied  as  theoretical  entities.  Traditional 
numerical  programs  are  hand  crafted  for  each  application.  The  traditional 
style  does  not  admit  such  explicit  decomposition  and  naming  of  methods, 
thus  losing  a  great  part  of  the  power  and  joy  of  programming. 

Appendix:  Convergence  rate  of  the  first-terms  sequence 

We  offer  here  the  mathematical  justification  for  claims  made  earlier  about 
the  rate  of  convergence  of  the  sequence  of  first  terms  returned  by  the  proce¬ 
dure  richardson-sequence. 

As  for  hypotheses,  we  suppose  that  R  is  a  function  analytic  in  a  disk  that 
contains  h  as  an  interior  point;  thus  we  have  an  expansion 

OO 

R(h)  =  A  +  ^Eih”, 

i=  1 

where  {pi,P2,-  •  •}  is  an  increasing  sequence  of  positive  integers.  We  assume 
the  pi  are  chosen  so  that  no  E{  is  0.  Let  us  identify  R  with  /?W;  the  sequence 
of  first  terms  is  given  by 

where  for  all  n  >  0, 

1  ’  2Pn  —  1 

This  is  the  operation  concocted  to  remove  the  dominant  error  term  at 
each  stage;  thus  we  know 


Since  it  is  possible  that  the  coefficients  grow  large  as  n  — ►  oo,  we  cannot 
immediately  conclude  that  Sn  converges  to  A  or  even  converges  at  all.  We 
settle  this  question  by  a  straight-forward  computation. 

To  begin  with,  we  have 


S2  —  A  + 


1 

2P1  -  1 


£«•( 


«=2 


1 

2w-pi 


-  1)A* 


and 


1  00  |  2 

'^3  =  A  +  (2pi  -  1)(2ps  -  1)  ^^'^Pi-Pi  _  ^^p.-w  _  l^hP  ' 
the  general  case  is  seen  to  be 

The  term  appearing  in  braces  is  less  than  1  in  magnitude;  this  gives  us 
the  estimate 

IS*.  -  ■*!  S  (ft  srrr)  £  IWI- 

\'=1  ^  1  /  i=n+l 

Since  R  is  analytic  at  h,  the  summation  part 

OO 

/<»=  £  \E,h”‘\ 

i=n-f  1 

converges  monotonically  to  0  as  n  — »  00.  The  product  appearing  in  paren¬ 
theses  is  estimated  as  follows. 

Let  cr„  =  E"=i  Pi\  then 

n— — = — nf — 

fA[  2V>  -  1  2°"  fJi  \  1  -  2-p*  ) 

Using  the  inequality,  valid  for  x  €  [0,  |], 


h 


we  see  that  the  product  above  is  dominated,  in  magnitude,  by 


where  K  =  n^oU  +  2~‘)  is  an  absolute  constant. 

Thus  we  have  shown  that  the  absolute  error  with  which  5n+1  approxi¬ 
mates  the  limit  A  is  less  than 

l<h » 

2a" 

In  the  cases  cited  in  our  discussion,  {p,}  is  an  arithmetic  progression  {p  + 
(:  -  !)<?};  hence 

\ ^  /  /•  ,\\  Tl{Tl  11 

an  =  2^(p  +  (*  -  1)9)  =  np+ - - - q. 

i= 1  1 

This  justifies  our  earlier  claim. 
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